QTL2 Mediation

Brian Yandell

library(ggplot2)

Specific analysis setup

These set up the phenotype and its peak location for further analysis.

chr_id <- "11"
scan_window <- c(96,98)
pheno_name <- "LD_light_pct"
peak_Mbp <- mean(scan_window)
window_Mbp <- diff(scan_window) / 2

Setup

local <- TRUE # consider local QTL (from mRNA traits)
other <- FALSE # consider other types of phenotypes if TRUE
qtls <- 1 # or 2 for number of drivers

The qtl2shiny package uses a data frame of project information (folder locations).

dirpath <- "../../qtl2shinyApp"
project_info <- data.frame(taxa = "CCmouse",
                           project = "Recla",
                           directory = "qtl2shinyData")

These encode a hierarchy:

datapath <- file.path(file.path(dirpath, project_info$directory))
taxapath <- file.path(datapath, project_info$taxa)
projectpath <- file.path(taxapath, project_info$project)

Information for all phenotypes: all covariates and the physical map pmap.

covar <- readRDS(file.path(projectpath, "covar.rds"))
str(covar)
## 'data.frame':    261 obs. of  6 variables:
##  $ sex       : chr  "male" "male" "male" "male" ...
##  $ Cohort    : chr  "Fall" "Fall" "Fall" "Fall" ...
##  $ Group     : chr  "1" "1" "1" "1" ...
##  $ Subgroup  : chr  "1B" "1A" "1A" "1A" ...
##  $ ngen      : chr  "4" "4" "4" "4" ...
##  $ coat_color: chr  "AG" "AG" "AG" "WH" ...
pmap <- readRDS(file.path(projectpath, "pmap.rds"))
sapply(pmap, range)
##               1          2          3          4          5          6
## [1,]   3.252796   3.164247   3.163885   3.427907   3.160082   3.278808
## [2,] 194.886567 181.812247 159.490311 156.138564 151.428385 149.460545
##               7          8          9         10         11         12
## [1,]   3.215308   3.235539   3.582987   3.180584   3.355458   3.562561
## [2,] 145.309341 128.662187 123.970379 130.550456 121.946117 119.342990
##              13         14         15       16        17        18        19
## [1,]   4.466555   3.262695   3.269262  3.16841  3.264958  3.235759  3.286882
## [2,] 119.392941 124.681199 103.176510 97.89090 94.717592 90.441479 61.167435
##               X
## [1,]   5.503248
## [2,] 169.979596

Now subset on region of interest:

pmap <- pmap[[chr_id]]

These are used to find objects that are not immediately loaded, including genotype probabilities and gene and variant information. See next section on query routines.

The phenotype data are extracted here.

pheno_data <- readRDS(file.path(projectpath, "pheno_data.rds"))
dim(pheno_data)
## [1] 261  26
colnames(pheno_data)
##  [1] "OF_distance_first4"        "OF_distance"              
##  [3] "OF_corner_pct"             "OF_periphery_pct"         
##  [5] "OF_immobile_pct"           "OF_center_pct"            
##  [7] "LD_distance_light"         "LD_transitions"           
##  [9] "LD_light_first4"           "LD_light_pct"             
## [11] "VC_top_distance_first4"    "VC_top_time_first4"       
## [13] "VC_top_distance"           "VC_top_time_pct"          
## [15] "VC_top_velocity"           "VC_bottom_distance_first4"
## [17] "VC_bottom_time_first4"     "VC_bottom_distance"       
## [19] "VC_bottom_time_pct"        "VC_bottom_velocity"       
## [21] "VC_bottom_transitions"     "HP_latency"               
## [23] "TS_frequency_climbing"     "TS_time_immobile"         
## [25] "TS_latency_immobile"       "bw"

The peaks and analyses_tbl have information on phenotypes. The peaks has location of all previously found peaks, used to find what maps in a region. The analyses_tbl has details of transformations and covariates. Phenotypes are ordered into groups and types, which is useful for searching and for determining sets of covariates and transformations when there are many phenotypes to consider.

peaks <- readRDS(file.path(projectpath, "peaks.rds"))
str(peaks)
## 'data.frame':    450 obs. of  8 variables:
##  $ pheno      : chr  "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
##  $ chr        : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ pos        : num  84.1 115.9 26.5 21.7 128.5 ...
##  $ lod        : num  5.25 5.67 6.55 3.65 5.03 ...
##  $ longname   : chr  "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
##  $ output     : chr  "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
##  $ pheno_group: chr  "group" "group" "group" "group" ...
##  $ pheno_type : chr  "type" "type" "type" "type" ...
analyses_tbl <- readRDS(file.path(projectpath, "analyses.rds"))
str(analyses_tbl)
## 'data.frame':    26 obs. of  15 variables:
##  $ pheno      : chr  "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
##  $ longname   : chr  "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
##  $ output     : chr  "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
##  $ pheno_group: chr  "group" "group" "group" "group" ...
##  $ pheno_type : chr  "type" "type" "type" "type" ...
##  $ model      : chr  "normal" "normal" "normal" "normal" ...
##  $ transf     : chr  "identity" "identity" "identity" "identity" ...
##  $ offset     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ winsorize  : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ sex        : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ Cohort     : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ Group      : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ Subgroup   : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ ngen       : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ coat_color : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

Query routines

The query routines embed the location and database information for genes, variants, probabilities and mRNA data, which are generally very large and only needed for a subset. All have first three arguments (chr, start, end)

They are used in read_project and read_query_rds, but we make their use explicit here for encoding locations in query functions.

query_genes <- qtl2::create_gene_query_func(file.path(taxapath, "mouse_genes.sqlite"))
query_variants <- qtl2::create_variant_query_func(file.path(taxapath, "cc_variants.sqlite"))
args("query_genes")
## function (chr, start, end) 
## NULL
str(query_genes(chr_id, scan_window[1], scan_window[2]))
## 'data.frame':    11201 obs. of  15 variables:
##  $ chr    : chr  "11" "11" "11" "11" ...
##  $ source : chr  "VEGA" "VEGA" "ENSEMBL" "ENSEMBL" ...
##  $ type   : chr  "gene" "mRNA" "gene" "mRNA" ...
##  $ start  : num  96 96 96 96 96 ...
##  $ stop   : num  96 96 96 96 96 ...
##  $ score  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ strand : chr  "-" "-" "-" "-" ...
##  $ phase  : num  NA NA NA NA NA NA NA NA 2 NA ...
##  $ ID     : chr  "VEGA:OTTMUSG00000001937" "VEGA:OTTMUST00000003877" "ENSEMBL:ENSMUSG00000013415" "ENSEMBL:ENSMUST00000013559" ...
##  $ Name   : chr  NA NA NA NA ...
##  $ Parent : chr  NA "VEGA:OTTMUSG00000001937" NA "ENSEMBL:ENSMUSG00000013415" ...
##  $ Dbxref : chr  "MGI:MGI:1890357" "MGI:MGI:1890357,VEGA:OTTMUSG00000001937" "MGI:MGI:1890357" "MGI:MGI:1890357,ENSEMBL:ENSMUSG00000013415" ...
##  $ mgiName: chr  NA NA NA NA ...
##  $ bioType: chr  "KNOWN_protein_coding" NA "protein_coding" NA ...
##  $ Alias  : chr  NA NA NA NA ...
args("query_variants")
## function (chr, start, end) 
## NULL
str(query_variants(chr_id, scan_window[1], scan_window[2]))
## 'data.frame':    25588 obs. of  16 variables:
##  $ snp_id      : chr  "rs218229099" "rs241313271" "rs27071476" "rs220108253" ...
##  $ chr         : chr  "11" "11" "11" "11" ...
##  $ pos         : num  96 96 96 96 96 ...
##  $ alleles     : chr  "A|T" "G|T" "A|G" "G|A" ...
##  $ sdp         : int  225 225 128 225 97 32 128 97 128 64 ...
##  $ ensembl_gene: chr  "ENSMUSG00000013415" "ENSMUSG00000013415" "ENSMUSG00000013415" "ENSMUSG00000013415" ...
##  $ consequence : chr  "intron_variant" "intron_variant" "intron_variant" "intron_variant" ...
##  $ A_J         : num  2 2 1 2 2 1 1 2 1 1 ...
##  $ C57BL_6J    : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ 129S1_SvImJ : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ NOD_ShiLtJ  : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ NZO_HlLtJ   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ CAST_EiJ    : num  2 2 1 2 2 2 1 2 1 1 ...
##  $ PWK_PhJ     : num  2 2 1 2 2 1 1 2 1 2 ...
##  $ WSB_EiJ     : num  2 2 2 2 1 1 2 1 2 1 ...
##  $ type        : chr  "snp" "snp" "snp" "snp" ...

The query_probs assumes a probdir (default genoprob), which contains the FST method files probs_fstindex.rds (for genotype probabilities) and aprobs_fstindex.rds (for allele probabilities). These are FST master files, which point to individual FST databases by chromosome, such as aprobs_1.fst and probs_1.fst.

query_probs <- qtl2pattern::create_probs_query_func(projectpath)
args("query_probs")
## function (chr = NULL, start = NULL, end = NULL, allele = TRUE, 
##     method = method_val, probdir = probdir_val) 
## NULL
get("method_val", env = environment(query_probs))
## [1] "fst"
get("probdir_val", env = environment(query_probs))
## [1] "genoprob"
probs_obj <- query_probs(chr_id, scan_window[1], scan_window[2], allele = TRUE)
dim(probs_obj$probs[[1]])
## [1] 261   8   9
range(probs_obj$map)
## [1] 96.04762 97.95779
dim(query_probs(chr_id, scan_window[1], scan_window[2], allele = FALSE)$probs[[1]])
## [1] 261  36   9

The query_mrna assumes a mrnadir (default RNAseq), which should contain the following files:

annot.mrna.rds
annot.samples.rds
expr.mrna.fst
peaks.mrna.fst
query_mrna <- qtl2mediate::create_mrna_query_func(projectpath)

The function has additional arguments local (use only mRNA whose genes are in the region if TRUE) and qtl (use only mRNA whose QTL peaks are in the region if TRUE).

args("query_mrna")
## function (chr = NULL, start = NULL, end = NULL, local = TRUE, 
##     qtl = FALSE, mrnadir = mrnadir_val) 
## NULL
get("mrnadir_val", env = environment(query_mrna))
## [1] "RNAseq"

There are no mRNA data for the Recla dataset, so this will not be used below. The following will create an error for these data.

str(query_mrna(chr_id, scan_window[1], scan_window[2]))

Phenotype Mediation

Phenotype mediation uses as mediators other phenotypes in the region defined by chr_id and scan_window. This is all bundled in the routine comediator_region, using project_info and analyses_tbl settings. It then calls pheno_region with peaks and qtls.

comed_ls <- qtl2mediate::comediator_region(pheno_name, chr_id, scan_window, 
                      covar, analyses_tbl, peaks, 
                      qtls, pmap, pheno_data)
med_ls <- qtl2mediate:::comediator_type(comed_ls, peaks, pheno_name, other)

Expression Mediation

Need to either have query_mrna set up ahead of time (in ‘RDS’ format) or create here.

expr_ls <- get_expr_region(chr_id, scan_window[1], scan_window[2], covar, pmap, 
                           drivers = qtls,
                           query_mrna = query_mrna)

Mediation Test

Mediation test is performed in shinyMediate. Input includes med_ls (from other phenotypes) or expr_ls (from mRNA expression data). Genetic cross data has kinship information, which is important to adjust for genetic correlation among individuals.

kinship <- readRDS(file.path(projectpath, "kinship.rds"))[[chr_id]]

The position for target mediation is by default the center of the map in this region. However, it might be improved by looking at peak for pheno_name.

pos_Mbp <- mean(range(probs_obj$map))
peak_mar <- qtl2::find_marker(probs_obj$map, chr_id, pos_Mbp)
geno_max <- subset(probs_obj$probs, chr = chr_id, mar = peak_mar)[[1]][,,1]
cov_tar <- qtl2mediate:::get_covar(
  covar, 
  dplyr::filter(analyses_tbl, pheno == pheno_name))
driver_med <- probs_obj$probs[[chr_id]]
phe_mx <- pheno_data[, pheno_name, drop = FALSE]
med_scan <- intermediate::mediation_scan(
  target = phe_mx,
  mediator = med_ls[[1]],
  driver = geno_max,
  annotation = med_ls[[2]],
  covar = cov_tar,
  kinship = kinship,
  fitFunction =  qtl2mediate::fitQtl2)
summary(med_scan, minimal = TRUE)
##                        chr      pos        LR
## LD_distance_light       11 97.02702  3.188115
## LD_transitions          11 96.68386  8.798000
## VC_bottom_distance      11 97.02702  9.884018
## VC_top_distance_first4  11 96.17147 10.321125
## OF_distance             11 97.02702 11.709825
## VC_bottom_transitions   11 97.29928 12.331006
## OF_distance_first4      11 97.02702 12.596732
ggplot2::autoplot(med_scan, size = 3)

med_test <- qtl2mediate::mediation_test_qtl2(
  target = phe_mx,
  mediator = med_ls[[1]],
  annotation = med_ls[[2]],
  covar_tar = cov_tar,
  covar_med = med_ls$covar,
  genoprobs = probs_obj$probs,
  map = probs_obj$map,
  chr = chr_id,
  pos = pos_Mbp,
  kinship = kinship)
summary(med_test)
## # A tibble: 7 × 30
##   id      chr     pos triad pvalue alt_t…¹ undec…² LR_me…³ LR_me…⁴ LR_CMST    df
##   <chr>   <fct> <dbl> <fct>  <dbl> <fct>     <dbl>   <dbl>   <dbl>   <dbl> <dbl>
## 1 LD_dis… 11     97.0 caus…      0 reacti…       0    3.19    15.1   141.      8
## 2 LD_tra… 11     96.7 caus…      0 reacti…       0    8.8     13      31.5     8
## 3 OF_dis… 11     97.0 caus…      0 reacti…       0   11.7     10.6    21.8     8
## 4 OF_dis… 11     97.0 caus…      0 reacti…       0   12.6     11.5    20.0     8
## 5 VC_bot… 11     97.0 caus…      0 reacti…       0   11.1     11.7    28.6     8
## 6 VC_bot… 11     97.3 caus…      0 reacti…       0   13.5      8.5    20.8     8
## 7 VC_top… 11     96.2 caus…      0 reacti…       0   11.5     11.1    18.8     8
## # … with 19 more variables: IC <dbl>, lod <dbl>, longname <chr>, output <chr>,
## #   pheno_group <chr>, biotype <chr>, model <chr>, transf <chr>, offset <dbl>,
## #   winsorize <lgl>, sex <lgl>, Cohort <lgl>, Group <lgl>, Subgroup <lgl>,
## #   ngen <lgl>, coat_color <lgl>, qtl_ct <int>, info <chr>, local <lgl>, and
## #   abbreviated variable names ¹​alt_test, ²​undecided, ³​LR_mediation,
## #   ⁴​LR_mediator
ggplot2::autoplot(med_test, size = 3, alpha = 1)

Mediation Triad

Based on shinyTriad in qtl2shiny. Basic idea is to find traits with strongest additive allele and/or SNP signal in the region based on scans. Then uses these to identify the best SDP(s) for further investigation with triads.

snpinfo <- query_variants(chr_id, scan_window[1], scan_window[2])
snpprobs <- qtl2mediate:::get_snpprobs(chr_id, peak_Mbp, window_Mbp,
                   pheno_name, 
                   probs_obj$probs,
                   probs_obj$map,
                   snpinfo)
snp_scan_obj <- qtl2mediate:::scan1covar(pheno_data, cov_tar,
                                        snpprobs$snpprobs, kinship,
                                        analyses_tbl,
                                        sex_type = "all")
qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
                         lodcolumn = pheno_name,
                         patterns = "hilit", drop_hilit = 3,
                         haplos = snpprobs$haplos,
                         main = pheno_name)

Want to use qtl2::top_snps instead of qtl2pattern::top_snps_pattern, but then have to do the work. Put this in triad_helper() in ’medation_triad_qtl2.R`.

(tmp_snps_1 <- 
  dplyr::arrange(
    qtl2::top_snps(snp_scan_obj, 
                   snpprobs$snpinfo,
                   lodcolumn = pheno_name,
                   show_all_snps = FALSE),
    dplyr::desc(.data$lod)))
##        snp_id chr      pos alleles sdp ensembl_gene        consequence A_J
## 1 rs257161856  11 97.00586     G|A 116              intergenic_variant   1
## 2  rs27023987  11 97.00335   C|G/T 139              intergenic_variant   3
##   C57BL_6J 129S1_SvImJ NOD_ShiLtJ NZO_HlLtJ CAST_EiJ PWK_PhJ WSB_EiJ type index
## 1        1           2          1         2        2       2       1  snp 12491
## 2        2           1          2         1        1       1       2  snp 12464
##   interval on_map      lod
## 1        3  FALSE 4.884181
## 2        3  FALSE 4.884181
qtl2mediate:::sdp_to_pattern(tmp_snps_1$sdp[1:2], haplos = snpprobs$haplos)
## [1] "ABDH:CEFG" "ABDH:CEFG"

Not quite what we want here. Want top SDP and top peaks; 1 or more.

snp_id_sdp <- dplyr::filter(snpprobs$snpinfo, sdp %in% tmp_snps_1$sdp[1:2])$snp_id
m <- match(snp_id_sdp, rownames(snp_scan_obj), nomatch = 0)
tmp <- as.data.frame(t(snp_scan_obj[m,]))
dplyr::arrange(tmp, dplyr::desc(.data[[names(tmp)[1]]]))
##                            rs27023987 rs257161856
## LD_light_pct              4.884180816 4.884180816
## LD_distance_light         4.290708137 4.290708137
## LD_transitions            3.536163536 3.536163536
## VC_top_distance_first4    3.533595886 3.533595886
## VC_top_distance           2.592797951 2.592797951
## OF_immobile_pct           1.932726913 1.932726913
## VC_top_velocity           1.486847315 1.486847315
## OF_distance               1.042637899 1.042637899
## VC_bottom_distance        1.006557071 1.006557071
## TS_frequency_climbing     0.952787113 0.952787113
## HP_latency                0.717622729 0.717622729
## OF_corner_pct             0.702872663 0.702872663
## VC_bottom_transitions     0.584050714 0.584050714
## OF_distance_first4        0.550941369 0.550941369
## VC_bottom_velocity        0.503672483 0.503672483
## TS_time_immobile          0.441407270 0.441407270
## LD_light_first4           0.406734913 0.406734913
## OF_center_pct             0.392502502 0.392502502
## bw                        0.359302531 0.359302531
## VC_bottom_distance_first4 0.292490619 0.292490619
## OF_periphery_pct          0.224770141 0.224770141
## VC_top_time_pct           0.147379093 0.147379093
## VC_bottom_time_pct        0.143988264 0.143988264
## TS_latency_immobile       0.143689933 0.143689933
## VC_bottom_time_first4     0.008482029 0.008482029
## VC_top_time_first4        0.003921065 0.003921065
qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
                         lodcolumn = pheno_name,
                         patterns = "hilit", drop_hilit = 3,
                  haplos = snpprobs$haplos)

qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
                         lodcolumn = match(c("LD_light_pct",
                                             "LD_transitions", 
                                             "LD_distance_light"),
                                           colnames(snp_scan_obj)),
                         #seq_len(ncol(snp_scan_obj))[1:2],
                         patterns = "hilit", drop_hilit = 3,
                         facet = "pheno",
                  haplos = snpprobs$haplos)

top_snps_tbl <- qtl2pattern::top_snps_pattern(
  snp_scan_obj,
  snpprobs$snpinfo)
summary(top_snps_tbl)
## # A tibble: 36 × 8
##    pheno                  min_pos max_pos max_lod min_lod   sdp pattern   snp_id
##    <chr>                    <dbl>   <dbl>   <dbl>   <dbl> <int> <chr>     <chr> 
##  1 LD_light_pct              97.0    97.0    4.88    4.88   116 ABDH:CEFG rs257…
##  2 LD_light_pct              97.0    97.0    4.88    4.88   139 ABDH:CEFG rs270…
##  3 LD_transitions            96.7    96.8    4.46    3.42    96 ABCDEH:FG 1819 …
##  4 VC_top_distance_first4    96.7    96.7    4.36    3.82   138 BDH:ACEFG 19 SN…
##  5 LD_distance_light         97.0    97.0    4.29    4.29   116 ABDH:CEFG rs257…
##  6 LD_distance_light         97.0    97.0    4.29    4.29   139 ABDH:CEFG rs270…
##  7 LD_distance_light         97.0    97.0    3.86    3.40   101 BDEH:ACFG 8 SNPs
##  8 VC_top_distance_first4    96.7    96.7    3.82    3.82   117 BDH:ACEFG rs249…
##  9 VC_top_distance_first4    96.6    96.7    3.74    3.74   101 BDEH:ACFG 6 SNPs
## 10 LD_transitions            97.0    97.0    3.54    3.54   116 ABDH:CEFG rs257…
## # … with 26 more rows
patterns <- summary(top_snps_tbl)
haplos <- LETTERS[1:8]

triad <- unique(med_test$best$triad)

# Pick top trait with causal call.
med_name <- dplyr::filter(med_test$best, .data$triad == "causal")[["id"]][1]
pattern <- dplyr::filter(patterns, .data$pheno == med_name)$pattern[1]

sdps <- unique(dplyr::filter(patterns, .data$pheno == med_name)$sdp)

sdp <- sdps[qtl2mediate:::sdp_to_pattern(sdps, haplos) == pattern]
id <- med_ls[[2]]$id[med_ls[[2]][["id"]] == med_name]
med_triad <- intermediate::mediation_triad(target = phe_mx,
                              mediator = med_ls[[1]][, id, drop = FALSE],
                              driver = geno_max, 
                              covar_tar = cov_tar,
                              covar_med = med_ls$covar,
                              kinship = kinship,
                              sdp = sdp, allele = TRUE,
                              fitFunction = qtl2mediate::fitQtl2)

Mediation Triad with Covar routine

Look carefully at code in last section, and at triad_sdp. Is the goal to find the sdp for the best triad across mediators or is it to just look at preselected mediator?

sdp <- qtl2mediate::triad_sdp(chr_id, pos_Mbp,
                 scan_window,
                 pheno_name, pheno_data = pheno_data,
                 covar_tar = cov_tar,
                 probs_obj$probs, probs_obj$map,
                 analyses_tbl, 
                 kinship = kinship)

Need to customize to particular mediator chosen.

med_triad <- qtl2mediate::mediation_triad_qtl2(
  target = phe_mx,
  mediator = med_ls[[1]][, id, drop = FALSE],
  annotation = med_ls[[2]],
  covar_tar = cov_tar,
  covar_med = med_ls$covar,
  genoprobs = probs_obj$probs,
  map = probs_obj$map,
  chr = chr_id,
  pos = pos_Mbp,
  kinship = kinship,
  sdp)
summary(med_triad)
## Model fit with driver columns
## 
## |term             |    estimate| std.error|  statistic|   p.value|
## |:----------------|-----------:|---------:|----------:|---------:|
## |mediator         |   0.0062485| 0.0026272|  2.3784080| 0.0181488|
## |SexFemale        |  35.5354522| 5.6471836|  6.2925973| 0.0000000|
## |SexMale          |  33.7103908| 6.1341315|  5.4955442| 0.0000001|
## |A                |   2.0962236| 6.1885618|  0.3387255| 0.7351043|
## |B                |   0.7602089| 6.3443328|  0.1198249| 0.9047193|
## |C                | -15.6357536| 6.1863361| -2.5274659| 0.0121137|
## |D                |  -1.4510211| 5.9428114| -0.2441641| 0.8073064|
## |E                |  -4.7513767| 5.9949686| -0.7925607| 0.4287944|
## |F                | -14.8780078| 5.4475045| -2.7311603| 0.0067662|
## |G                | -14.7548869| 5.6738148| -2.6005232| 0.0098693|
## |H                |   0.0000000| 0.0000000|  0.0000000| 0.0000000|
## |mediator:SexMale |   0.0028075| 0.0038471|  0.7297690| 0.4662229|
## 
## Model fit with driver groups
## 
## |term             |    estimate|  std.error|  statistic|   p.value|
## |:----------------|-----------:|----------:|----------:|---------:|
## |mediator         |   0.0054653|  0.0028453|  1.9208470| 0.0560495|
## |SexFemale        |  42.7192755|  8.5793636|  4.9793059| 0.0000013|
## |SexMale          |  40.3536882|  8.6381221|  4.6715811| 0.0000052|
## |groupAB          |   0.8476695| 10.6113354|  0.0798834| 0.9364029|
## |groupAC          | -19.1925195|  8.8554365| -2.1673149| 0.0312895|
## |groupAD          | -21.1649362| 10.6412234| -1.9889570| 0.0479506|
## |groupAE          |  -6.5846890| 16.9331714| -0.3888633| 0.6977551|
## |groupAF          | -19.0213010|  8.8830746| -2.1412970| 0.0333546|
## |groupAG          |  -4.6151459|  8.5793530| -0.5379364| 0.5911672|
## |groupAH          |  -5.1486070|  8.7570079| -0.5879413| 0.5571774|
## |groupBB          |   2.9607322| 16.7707948|  0.1765410| 0.8600321|
## |groupBC          |  -5.1323361|  9.3957988| -0.5462373| 0.5854588|
## |groupBD          |  -3.7704978| 16.9150829| -0.2229074| 0.8238153|
## |groupBE          | -16.7587658| 10.5788947| -1.5841698| 0.1145980|
## |groupBF          | -14.7930949|  8.2749856| -1.7876883| 0.0752090|
## |groupBG          | -16.7602294|  9.0182555| -1.8584780| 0.0644429|
## |groupBH          |  -5.2897217|  9.7078816| -0.5448894| 0.5863840|
## |groupCC          |   7.2480651| 16.7691641|  0.4322258| 0.6660026|
## |groupCD          | -11.8523945| 10.0588926| -1.1783001| 0.2399557|
## |groupCE          | -15.0919469|  8.8576031| -1.7038409| 0.0898292|
## |groupCF          | -29.5018343|  9.2084014| -3.2037954| 0.0015584|
## |groupCG          | -28.0035295|  8.9753249| -3.1200575| 0.0020514|
## |groupCH          |  -7.3006264|  9.4034348| -0.7763787| 0.4383631|
## |groupDD          |  -7.4664462|  9.6836985| -0.7710325| 0.4415187|
## |groupDE          |  -2.9445950|  8.6502080| -0.3404074| 0.7338760|
## |groupDF          | -17.6802963|  8.5688101| -2.0633315| 0.0402602|
## |groupDG          | -22.6306973| 10.0363535| -2.2548725| 0.0251303|
## |groupDH          |  -2.1399215|  9.3840554| -0.2280380| 0.8198294|
## |groupEE          | -11.2954916| 11.4532466| -0.9862262| 0.3251105|
## |groupEF          | -19.4211281|  9.2044763| -2.1099656| 0.0359964|
## |groupEG          | -19.2922730|  9.2171838| -2.0930767| 0.0374936|
## |groupEH          | -11.5139640| 16.7588549| -0.6870376| 0.4927859|
## |groupFF          | -12.8569264|  9.1678473| -1.4023932| 0.1622138|
## |groupFG          | -15.7006508|  8.2758211| -1.8971714| 0.0591210|
## |groupFH          | -15.0044632|  8.4325119| -1.7793587| 0.0765678|
## |groupGG          | -29.5473583| 11.4551739| -2.5793898| 0.0105518|
## |groupGH          |  -4.4104300| 10.6128505| -0.4155745| 0.6781283|
## |groupHH          | -23.3427702| 11.6163964| -2.0094674| 0.0457148|
## |mediator:SexMale |   0.0035596|  0.0040655|  0.8755636| 0.3822266|
intermediate::ggplot_mediation_triad(med_triad)

intermediate::ggplot_mediation_triad(med_triad, fitlines = "parallel")

plotly::ggplotly(intermediate::ggplot_mediation_triad(med_triad, fitlines = "parallel"))